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Gravitational waves from the sound of a first order phase transition 
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We report on the first 3-dimensional numerical simulations of first-order phase transitions in the 
early universe to include the cosmic fluid as well as the scalar field order parameter. We calculate the 
gravitational wave (GW) spectrum resulting from the nucleation, expansion and collision of bubbles 
of the low-temperature phase, paying particular attention to those sourced by the fluid. We find 
that the fluid continues to be a source of GWs long after the bubbles have merged, a new effect not 
taken into account in previous modelling of the GW source based on the envelope approximation. 
The kinetic energy of the fluid is in the form of compression waves: the main source of the GWs 
after a phase transition is therefore the sound the bubbles make. 
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In a hot Big Bang there were phase transitions in the 
early Universe [l|, |2| , which may well have been of first 
order; one major consequence of such a transition would 
be the generation of gravitational waves [3|-l8| . The elec- 
troweak transition in the Standard Model is known to be 
a cross-over |9l4l l| but it may be first order in minimal 
extensions of the Standard Model |12l4l7j. It is there- 
fore essential to properly characterise the expected power 
spectrum from first-order phase transitions. 

First order phase transitions proceed by the nucleation, 
growth, and merger of bubbles of the low temperature 
phase [sl, ll8M25| . The collision of the bubbles is a violent 
process, and both the scalar order parameter and the 
fluid of light particles generate gravitational waves. 

Numerical studies have been carried out of the be- 
haviour of bubbles in such a phase transition using spher- 
ically symmetric (l + l)-dimensional simulations ^23l i24j] . 
The calculation of the gravitational wave spectrum has 
been refined in the intervening years using the semi- 
analytic envelope approximation [5|, tzl, [sl, l26l428| . Fully 
3D simulations of the scalar field only have been carried 
out [29|, qualitatively supporting the envelope approxi- 
mation, and pointing out important gravitational wave 
production from the scalar field after the bubble merger. 

In a hot phase transition, the fluid plays an important 
role, firstly as a brake on the scalar field, and secondly 
as a source of gravitational waves itself. The fluid has 
generall y b een assumed to be incompressible and turbu- 
lent 



rall y b e( 
30M33|. 



An important question for the gravitational 
wave power spectrum is the validity of this modelling, 
which generally borrows from the Kolmogorov theory of 
non-relativist ic driven incompressible turbulence. 

In this Letter we report on the first fully three dimen- 
sional simulation of bubble nucleation involving a cou- 
pled field- fluid system. We make use of these simulations 
to calculate the power spectrum of gravitational radia- 
tion from a first-order phase transition. We find that the 
compression waves in the fluid - sound waves - continue 



to be an important source of gravitational waves for up 
to a Hubble time after the bubble merger has completed. 
This boosts the signal by the ratio of the Hubble time to 
the transition time, which can be orders of magnitude. 

The system describing the matter in the early universe 
consists of a relativistic fluid coupled to a scalar field, 
which acquires an effective potential 
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The rest-frame pressure p and energy density e are 



(1) 



e = 3aT' + V{<P,T)-T— (2) 

p = aT^-V{4>,T) (3) 

with a = (7r^/90)^, and g the effective number of rela- 
tivistic degrees of freedom contributing to the pressure 
at temperature T. The stress-energy tensor for a scalar 
field (j) and an ideal relativistic fluid U^ is 

+ [e+p]/7^/7^+^^> (4) 

where the metric convention is (— + ++)• The scalar 
field potential is included in the definition of p. We split 
dfjT^^ = (nonuniquely) into field and fluid parts with 
a dissipative term permitting transfer of energy between 
the scalar field and the fluid (5^ = 7^/7^ S^05^^ l23i]. 

Given these expressions, the equations of motion can 
be derived. For the field we have 
dV 
- ^ + ^'^ - ^ = ^^(^ + ^'^^^) (^) 

0(j) 

where W is the relativistic 7-factor and V^ is the fluid 
3- velocity, U^ = WV^ . For the fluid energy density E = 
We, contracting [d/j^T^^]^^^^ with Uu yields 

dV 

E + d^{EV') -^p[W + d^{WV')] - -Q-W{(^ + V'd^cfy 




FIG. 1: Snapshot of scalar field at t = 400Te"\ for ry = 0.2 
with the 'weak' transition parameters. Regions with < 1 
are transparent. See Fig. [2] for the energy density on a spatial 
slice at this timestep. 




FIG. 2: Slices of fluid energy density E/T^ at t = 400Te~\ 
t = 800 Te"^ and t = 1200 T^"^ respectively, for the ry = 0.2 
simulation. The slices correspond roughly to the end of the 
nucleation phase, the end the end of the initial coalescence 
phase and the end of the simulation. 



The equations of motion for the fluid momentum density 

Zi = W{e^p)Ui read 

dV 



The principal observable of interest to us is the power 
spectrum of gravitational radiation resulting from bubble 
collisions. To calculate this we use the procedure detailed 
in Ref. [3^ . We evolve the unprojected equation of mo- 
tion for the metric perturbations 



V^Uij = WnGirf^ 
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where rf^ = di 



and rlj 



f _ 



W'^ {e ^ p)ViVj . We recover 
the physical metric perturbations in momentum space 
hij{\<i) = \ij^imi)<i)uim{t,^), where Ai^-,/rn(k) is the pro- 
jector onto transverse, traceless symmetric rank 2 ten- 
sors. We are most interested in the metric perturbations 
sourced by the fluid, as the fluid shear stresses generally 
dominate over those of the scalar field, although it will 
be instructive to also consider both sources together. 



Having obtained the metric perturbations, the power 
spectrum per logarithmic frequency interval is 
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d\nk 
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We simulate the system on a cubic lattice of N^ = 1024^ 
points, neglecting cosmic expansion which is slow com- 
pared with the transition rate. The fluid is imple- 
mented as a three dimensional relativistic fluid [35| , with 
donor cell advection. The scalar and tensor fields are 
evolved using a leapfrog algorithm with a minimal sten- 
cil for the spatial Laplacian. Principally we used lat- 
tice spacing 6x = IT~^ and time step St = 0.1 T~^, 
where Tc is the critical temperature for the phase tran- 
sition. We have checked the lattice spacing dependence 
by carrying out single bubble self-collision simulations for 
L^ = 256^ T-^ at 6x = 0.5T-\ for which the value of 
Pgw at t = 2000T~^ increased by 10%, while the final 
total fluid kinetic energy increased by 7%. Simulating 
with St = 0.2 T-i resulted in changes of 0.3% and 0.2% 
to pGw and the kinetic energy respectively. 

Starting from a system completely in the symmetric 
phase, we model the phase transition by nucleating new 
bubbles according to the probability per unit volume and 
time P = Poexp(/3(t — to)). Having decided to nucleate 
in a suitable untouched region of the box, we insert a 
static bubble with a gaussian profile for the scalar field. 

We first studied a system with g = 34.25, 7 = 1/18, 
a = V^/72, To = Tc/V2 and A = 10/648; this allows 
comparison with previous (1 -h 1) and spherical studies 
of a coupled field- fluid system where the same parame- 
ter choices were used [23|. The transition in this case 
is relatively weak: in terms of ar, the ratio between 
the latent heat and the total thermal energy, we have 
o^Tn = 0.012 at the nucleation temperature Tn = 0.86 Tc. 
We therefore also performed simulations with 7 = 2/18 
and A = 5/648, for which ar^ = 0-10 at the nucleation 
temperature Tn = 0.8 Tc, which we refer to as an inter- 
mediate strength transition. 

For the nucleation process, we took /3 = 0.0125 Tc, 
Po = 0.01 and to = tend = 2000 T-^ The simulation 
volume allowed the nucleation of 100-300 bubbles, cor- 
responding to a mean spacing between bubbles of order 
100 T~^. Although the wall velocity is captured correctly, 
there is not enough time for the bubble walls to fully relax 
to their invariant scaling profile before colliding. Typi- 
cally, the peak velocity prior to collision is 20-30% below 
the scaling value for the deflagrations. 

For the weak transition we chose r] = 0.1, 0.2, 0.4 and 
0.6. The first gives a detonation with wall speed v^ c^ 
0.71, and the others weak deflagrations with v^ c:^ 0.44, 
0.24, and 0.15 respectively. The shock profiles are found 
in Figs. 2 and 3 of Ref. [22^]. The intermediate transition 
was simulated at 7^ = 0.4, for which the wall speed is 
v^ ~ 0.44, very close to the weak transition with r] = 0.2. 
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FIG. 3: Top: time series of [/(/, and Ui fTO)) . showing the 
progress of the phase transition; the curves for U(f) and Ui are 
individually identified for the 'intermediate' case. Bottom: 
time series of pGwR^^{e -\- p)~'^Uf (tf), showing the evolu- 
tion of the gravitational wave energy density relative to an 
estimate of the square of the final fluid shear stresses. 



Fig. [3] (top) shows the time evolution of two quantities 
U(h and /7f, defined so that 



ie + p)U^ 
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d^xT[^ and (e + p)[/f 
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(10) 

where e and p are the volume- averaged rest-frame energy 
density and pressure respectively. 

The squares of these quantities give an estimate of the 
size of the shear stresses of the field and the fluid relative 
to the background fluid enthalpy density, while /7f tends 
to the r.m.s. fluid velocity for /7f <C 1. We see that ^7^ 
grows and decays with the total surface area of the bub- 
bles of the new phase, while the mean fluid velocity grows 
with the volume of the bubbles, and then stays constant 
once the bubbles have merged. The slight decreasing 
trend in U^ for the intermediate transition is a numerical 
effect: first-order advection is known to be diffusive, but 
it is adequate over the period of these simulations. 

Fig. [3] (bottom) shows the GW energy density scaled 




FIG. 4: Gravitational wave power spectra during the phase 
transition, for the intermediate strength transition, from fluid 
only (black) and both fluid and field (grey). From bottom to 
top, the times are t = 600, 800, 1000, 1200 and UOOT~^. 
The red dotted line indicates the expected k~^ behaviour. 



by the final value of (e +p)^t/f and the average bubble 
size at collision R^^ estimated from Rl = L^/Ni^, where 
N\) is the number of bubbles in the simulation volume. 
An isolated expanding bubble emits no GWs, which ap- 
pear only when the bubbles collide. We see that the GW 
energy density rises linearly after the bubbles have fully 
merged, sourced by the perturbations in the fluid. The 
slopes are similar, as expected for a source proportional 

2 

to (e +p)/7f with scale i?* (see discussion below). Note 
that the GWs from detonations (77 = 0.1) behave simi- 
larly to those from deflagrations. 

In Fig. |4] we show the time development of the GW 
power spectrum as the intermediate strength phase tran- 
sition proceeds. We see that strong growth happens be- 
tween t = 600 T-^ and t = 1000 T"^ as the bubbles 
merge (see Fig. [3|). For t < 1000 T~^ there is evidence 
of the expected k~^ power spectrum, but it becomes less 
clear as the GW power continues to grow, sourced by 
the persistent fluid perturbations. At the shortest length 
scales, we see a ^w-dependent exponential fall-off. 

To establish the nature of these fluid perturbations, we 
show in Fig. [5] the time development of the longitudinal 
(compressive) and transverse (rotational) components of 
the fluid velocity power spectrum. At all times, it is 
clear that most of the fluid velocity is longitudinal, in- 
dicating that the perturbations are mostly compression 
waves. Turbulence generally develops in the transverse 
components, characterised by a power-law behaviour of 
the power spectrum. There is some evidence of an almost 
flat spectrum for k < 0.2, but its significance is unclear. 

We can now form a clearer picture of the fluid per- 
turbations and how the GWs are generated. Firstly, we 
note that the fluid perturbations are initially the form 
of a compression wave surrounding the growing bubble. 
The energy in this compression wave grows in proportion 




FIG. 5: Fluid velocity power spectra for the intermediate 
strength transition, separated into longitudinal (compres- 
sional) and transverse (rotational) components; shown in grey 
and black respectively. Times shown are the same as Fig. [H 



to the volume of the bubble R^^ and quickly outstrips the 
energy in the scalar field, which grows only as R^. The 
energy in the compression waves remains constant after 
the bubble merger is complete: as fluid velocities are gen- 
erally much less than one, their evolution is linear. 

The bubble collision generates gravitational waves, as 
predicted by the envelope approximation, and there is 
some evidence for the characteristic k~^ spectrum be- 
tween the bubble separation scale R^ and the high- 
frequency cut-off. The generation of GWs continues long 
after the merger is completed and the scalar field has re- 
laxed to its new equilibrium value. The GWs are sourced 
by the compression waves in the fluid. This source of 
gravitational radiation from a phase transition - sound - 
has not been appreciated before (except in Ref. [4]). 

We can estimate the resulting density of the gravita- 
tional waves, in terms of the lifetime of the sound waves 
Ts, considered below. The power spectrum is given by 



dpGwjk) _ 2Gk^ 
dliik 7T 



/ dtidt2Cos[k{ti -fe)]n2(A:,ti,t2), 

(11) 

where II^(/c, ti, ^2) is the unequal time correlator of the 
shear stress tensor [28|, 



(e + p)C/f 



We can estimate the ampli- 
and the length scale as R^. The 



tude as 

coherence time can be written {csk)~^fc{kR^). Hence 



U\k,ti,t2) ^ [{e + p)UfYRiU^{kR.,Csz/fc), (12) 

where z = k{ti — ^2). The density parameter flow = 
pGw/^ is, on using the Friedmann equation, 

^GW ^ ^{H,n){H,R,){l + wfut, (13) 
where H^ is the Hubble parameter at the transition, w = 



p/e ^ 1/3, and 

H^ = / d\nk{kR^f / dzcos{z)IV(kR^,c^z/fc). (14) 

In (p!3|) we see the origin of the R^ factor in the GW 
density, which must be present for dimensional reasons. 
The slope of the curves in Fig. (j3l bottom) is 211^ /tt, 
which we see takes the natural value 0(1), and is weakly 
dependent on the transition parameters. 
The envelope approximation gives |27j 
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where k is the efficiency with which latent heat is con- 
verted to kinetic energy. Gomparing to ([13]) and noting 

that /7f ~ n^a^^ R^ ~ '^w//^, we see that sound waves 
are a more important source by the factor r^/R^v^. 

An upper bound on Tg is the Hubble time, as the shear 
stresses decay faster than the background energy den- 
sity. The shear stresses also decay due to the viscosity r^g, 
which can be estimated as rj^ ^ T^/e^ ln(l/e), where e is 
the electromagnetic gauge coupling [38]. The lifetime of 
sound waves with characteristic wavelength R^ is there- 
fore r^ c^ Rle/% ^ e^ln(l/e)i?^Tc. Hence sound waves 
are the most important source of gravitational waves for 
bubbles satisfying 

R^H^ > v^{^/^Tc/mpe^) - 10-^^Vw(Tc/100 GeV). 

(16) 
This is generally satisfied except for weak transitions at 
very high temperatures. 

We have carried out three dimensional simulations of 
bubble nucleation, growth and collision during a first or- 
der phase transition in the early universe, calculating 
the resulting gravitational wave power spectrum. We 
have included a relativist ic fluid for the first time. We 
find that, irrespective of whether it proceeds by deflagra- 
tion or detonation, the transition creates a background 
of acoustic waves, which are a more important source of 
gravitational radiation by the ratio of the fluid damping 
time Tg to the duration of the phase transition /3~^. For 
a wide range of parameters, the fluid damping time is the 
Hubble time, and we conclude that the dominant source 
of the gravitational wave background from a first order 
phase transition is the sound it makes. 
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cility jointly funded by STFC and the Large Facilities 
Capital Fund of BIS). KR acknowledges support from 
the Academy of Finland project 1134018; MH and SH 
from the Science and Technology Facilities Council (grant 
number ST/J000477/1) 
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